Skip to main content

Full text of "BSTJ 48: 5. May-June 1969: The Chirp z-Transform Algorithm and Its Application. (Rabiner, Lawrence R.; Schafer, Ronald W. ; Rader, Charles M.)"

See other formats


The Chirp z-Transform Algorithm 
and Its Application 

S By LAWRENCE R. RABINER, RONALD W. SCHAFER, 
and CHARLES M. RADER* 

(Manuscript received November 21, 1968) 

We discuss a computational algorithm for numerically evaluating the z- 
transform of a sequence of N samples. This algorithm has been named the 
chirp z-transform algorithm. Using this algorithm one can efficiently 
evaluate the z-transform at M points in the z-plane which lie on circular 
or spiral contours beginning at any arbitrary point in the z-plane. The 
angular spacing of the points is an arbitrary constant; M and N are 
arbitrary integers. 

The algorithm is based on the fact that the values of the z-transform on a 
circular or spiral contour can be expressed as a discrete convolution. Thus 
one can use well-known high-speed convolution techniques to evaluate the 
transform efficiently. For M and N moderately large, the computation 
time is roughly proportional to (N + M) log 2 (N + M) as opposed to 
being proportional to N-M for direct evaluation of the z-transform at M 
points. 

Applications discussed include: enhancement of poles in spectral analysis, 
high resolution narrow-band frequency analysis, interpolation of band- 
limited waveforms, and the conversion of a base 2 fast Fourier transform 
program into an arbitrary radix fast Fourier transform program. 

I. INTRODUCTION 

In dealing with sampled data the 2-transform plays the role which 
is played by the Laplace transform in continuous time systems. One 
example of its application is spectrum analysis. The computation of 
sampled ^-transforms, which has been greatly facilitated by the fast 
Fourier transform algorithm, is further facilitated by the "chirp 
z-transform" algorithm described in this paper. 1 - 2 

* Mr. Rader is with Lincoln Laboratory, Massachusetts Institute of Tech- 
nology, Lexington, Massachusetts. Lincoln Laboratory is operated with support 
from the U. S. Air Force. 

1249 



1250 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

The z-transform of a sequence of numbers x n is denned as 

X(z) = 2 x n z~ n , (1) 



»- — CO 



a function of the complex variable z. In general both x n and X(z) 
could be complex. It is assumed that the sum on the right side of equa- 
tion (1) converges for at least some values of z. We restrict ourselves 
to the z-transform of sequences with only a finite number N of non- 
zero points. Therefore, we can rewrite equation (1) without loss of 
generality as 



AT-l 



X(z) = Z *** (2) 



i = 



where the sum in equation (2) converges for all z except z = 0. 

Equations (1) and (2) are similar to the defining expressions for 
the Laplace transform of a train of equally spaced impulses of magni- 
tudes x„. Let the spacing of the impulses be T and let the train of 
impulses be 

Y,x n b{t-nT). 
Then the Laplace transform is 

n 

which is the same as X (z) if we let 

z = e' T . (3) 

For sampled waveforms the relation between the original waveform 
and the train of impulses is well understood in terms of the phenomenon 
of aliasing. Thus the z-transform of the sequence of samples of a time 
waveform is representative of the Laplace transform of the original 
waveform in a way which is well understood. The Laplace transform of a 
train of impulses repeats its values taken in a horizontal strip of the 
s-plane of width 2ir/T in every other strip parallel to it. The z-transform 
maps each such strip into the entire z-plane or, conversely, the entire 
z-plane corresponds to any horizontal strip of the s-plane, for example, 
the region -co < a < a>, —r/T ^ a> < ir/T, where s = a -+- jo. 
In the same correspondence, the ;'co axis of the s-plane, along which 
we generally equate the Laplace transform with the Fourier transform, 
is the unit circle in the z-plane; the origin of the s-plane corresponds to 
z ■* 1. The interior of the z-plane unit circle corresponds to the left 



CHIRP Z-TRANSFORM 



1251 



half of the s-plane; the exterior corresponds to the right half plane. 
Straight lines in the s-plane correspond to circles or spirals in the z-plane. 
Figure 1 shows the correspondence of a contour in the s-plane to a 
contour in the z-plane. To evaluate the Laplace transform of the im- 
pulse train along the linear contour is to evaluate the z-transform of 
the sequence along the spiral contour. 

Values of the z-transform are usually computed along the path 
corresponding to the ;« axis, namely the unit circle. This gives the 
discrete equivalent of the Fourier transform and has many applications 
including the estimation of spectra, filtering, interpolation, and corre- 
lation. The applications of computing z-transforms off the unit circle 
are fewer, but one is presented in this paper, namely the enhancement of 
spectral resonances in systems for which one has some foreknowledge 
of the locations of poles and zeros. 

Just as we can only compute equation (2) for a finite set of samples, 
so we can only compute equation (2) at a finite number of points, 
say z fc : 



x k - x(z k ) = £ « 



f k \ 



(4) 



The special case which has received the most attention is the set of 
points equally spaced around the unit circle, 



->■■ - ( ' X I' U ' j^* 



k = 0, 1 



,N - I 



(5) 



z-PLANE 





S-PLANE 



Fig. 1 — The correspondence of a z-plane contour to an s-plane contour through 
the relation z = e' T . 



1252 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

for which 



X k = E^expl-jf^nft] , k = 0, 1, ••• ,N - 1. (6) 



n-(l 



Equation (6) is called the discrete Fourier transform. The reader may 
easily verify that, in equation (5), other values of k merely repeat 
the same N values of z k , which are the iVth roots of unity. The discrete 
Fourier transform has assumed considerable importance, partly be- 
cause of its nice properties but mainly because since 1965 it has be- 
come widely known that the computation of equation (6) can be 
achieved, not in the N 2 complex multiplications and additions called 
for by direct application of equation (6), but in something of the 
order of N log2 iV operations if N is a power of two, or N 2 m* opera- 
tions if the integers ?w< are the prime factors of N. Any algorithm 
which accomplishes this is called a fast Fourier transform. Much of 
the importance of the fast Fourier transform is that the discrete 
Fourier transform may be used as a stepping stone to computing 
lagged products such as convolutions, autocorrelations, and cross cor- 
relations more rapidly than before. 3 - 4 The discrete Fourier transform 
has some limitations which can be eliminated using the chirp 2-trans- 
form algorithm which we describe. We also investigate the computa- 
tion of the 2-transform on a more general contour, of the form 

z k = AW~ k , k = 0, 1, • • • , M - 1 (7a) 

where M is an arbitrary integer and both A and W are arbitrary 
complex numbers of the form 

A = A exp O'27r0 o ) (7b) 

and 

W = "Pr o exp0"27r^). (7c) 

(See Fig. 2.) The case A = 1, M = N, and W = exp (—j2ir/N) corre- 
sponds to the discrete Fourier transform. The general z-plane contour 
begins with the point z = A and, depending on the value of W, spirals 
in or out with respect to the origin. If W = 1, the contour is an arc of a 
circle. The angular spacing of the samples is 2ir(p . The equivalent 
s-plane contour begins with the point 

s = o + jw = ^ In A (8) 



CHIRP Z -TRANSFORM 



1253 



S-PLANE 




Acr=y,lnW 



CT = y/nA 



Fig. 2 — An illustration of the independent parameters of the chirp z-transform 
algorithm. The upper figure shows how the z-transform is evaluated on a spiral 
contour starting at the point z = A. The lower figure shows the corresponding 
straight line contour and independent parameters in the s-plane. 

and the general point on the s-plane contour is 

s k = s + 7c(Acr + ; Aco) = ^(ln A - k In W), 

k = 0, 1, ••• ,M - 1. (9) 

Since A and W are arbitrary complex numbers we see that the points 
s fc lie on an arbitrary straight line segment of arbitrary length and 
sampling density. Clearly the contour indicated in equation (7a) is 
not the most general contour, but it is considerably more general 
than that for which the discrete Fourier transform applies. In Fig. 2, 
an example of this more general contour is shown in both the z-plane 
and the s-plane. 

To compute the z-transform along this more general contour would 
seem to require NM multiplications and additions since the special 
symmetries of exp (j2irk/N) which are exploited in the derivation of 
the fast Fourier transform are absent in the more general case. How- 
ever, we shall see that by using the sequence TF nV2 in various roles 
we can apply the fast Fourier transform to the computation of the 
z-transform along the contour of equation (7a). Since for W = 1, the 
sequence W n ' /2 is a complex sinusoid of linearly increasing frequency, 
and since a similar waveform used in some radar systems has the 
picturesque name "chirp", we call the algorithm we are about to present 
the chirp z-transform. Since this transform permits computing the 
z-transform on a more general contour than the fast Fourier transform 



1254 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

permits, it is more flexible than the fast Fourier transform, although it 
is also considerably slower. The additional freedoms offered by the 
chirp 2- transform include: 

(i) the number of time samples does not have to equal the number 
of samples of the z-transform. 

(it) Neither M nor N need be a composite integer. 
(Hi) The angular spacing of the z k is arbitrary. 
(iv) The contour need not be a circle but can spiral in or out with 
respect to the origin. In addition, the point z„ is arbitrary, but this is 
also the case with the fast Fourier transform if the samples x n are multi- 
plied by z~ n before transforming. 

II. DERIVATION OF THE CHIRP Z-TRANSFORM 

Along the contour of equation (7a), equation (4) becomes 

X h - S x n A' n W nk , h - 0, 1, ■ •• , M - 1 (10) 

n-0 

which, at first appearance, seems to require NM complex multiplica- 
tions and additions, as we have already observed. But, let us use 
Bluestein's ingenious substitution 5 

nk = »' + *'-(*- »)' (11) 

for the exponent of W in equation (10). This produces an apparently 
more unwieldly expression 

x k = T l x n A- n w (n,/2) w ik ' /2) w-«- n) ' /2 , 

n-0 

k = 0, 1, ••• ,M - 1 (12) 

but in fact equation (12) can be thought of as a three step process 
consisting of: (i) forming a new sequence y n by weighting the x n ac- 
cording to the equation 

y n = x n A-'W" /a , n - 0, 1, • • • , N - I, (13) 

(u) convolving y„ with the sequence v n defined as 

Vn = w~ n ' /2 (14) 

to give a sequence g k 



g» 



= Ep.- fc = 0,l, ■•• ,M-1, (15) 



CHIRP Z -TRANSFORM 



1255 



and (Hi) multiplying g k by W*' /a to give X k 



X k = g k W 



k = 0, 1, • ■ • , M - 1. 



(16) 



The three step process is illustrated in Fig. 3. Steps i and Hi require 
N and M multiplications respectively; step ii is convolution which may 
be computed by the high speed technique disclosed by Stockham, based 
on the use of the fast Fourier transform. 3 Step ii is the major part of 
the computational effort and requires a time roughly proportional to 
(N + M) log (N + M). 

Bluestein used the substitution of equation (11) to convert a dis- 
crete Fourier transform to a convolution as in Fig. 3. The linear sys- 
tem to which the convolution is equivalent can be called a chirp filter 
which is, in fact, also sometimes used to resolve a spectrum. Blue- 
stein showed that for N a perfect square, the chirp filter could be 
synthesized recursively with N^ multipliers and the computation of a 
discrete Fourier transform could then be proportional to N 3/2 (see 
Ref. 5) . 

The flexibility and speed of the chirp 2-transform algorithm are 
related to the flexibility and speed of the method of high-speed con- 
volution using the fast Fourier transform. Recall that the product of 
the discrete Fourier transforms of two sequences is the discrete Fourier 
transform of the circular convolution of the two sequences; therefore, 
a circular convolution is computable as two discrete Fourier trans- 
forms, the multiplication of two arrays of complex numbers, and an 
inverse discrete Fourier transform, which can also be computed by 
the fast Fourier transform. Ordinary convolutions can be computed 
as circular convolutions by appending zeroes to the end of one or 
both sequences so that the correct numerical answers for the ordinary 
convolution can result from a circular convolution. 

We now summarize the details of the chirp z-transform algorithm 
on the assumption that an already existing fast Fourier transform 



n = o,i,...,N-i 



<*> 



yn 



,nZ/2 



HIGH SPEED 
CONVOLUTION 
,-r>2/2 



y n ®w~ 



9k 



x k 



k = o,i,2,...M-i 



W 



k 2 /z 



Fig. 3 — An illustration of the steps involved in computing values of the 
z-transform using the chirp z-transform algorithm. 



1256 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

program (or special purpose machine) is available to compute discrete 
and inverse discrete Fourier transforms. 

We begin with a waveform in the form of N samples x n and we 
seek M samples of X k where A and W have also been chosen: 

(i) We choose L, the smallest integer greater than or equal to N + 
M — 1 which is also compatible with our high speed fast Fourier 
transform program. For most users this will mean L is a power of two. 
Notice that while many fast Fourier transform programs will work 
for arbitrary L, they are not equally efficient for all L. At the very 
least, L should be highly composite. 

(ii) We form an L point sequence y n from x n by the equation 

\A- n W n ' /2 x n n = 0, 1, 2, ■ • • , N - 1 n ~ 

y n = \ • (17) 

lo n = N,N + 1, ••• ,L - 1 

(Hi) We compute the L point discrete Fourier transform of y n by 
the fast Fourier transform, calling it Y r , r = 0, 1, . . . , L — 1. 
(iv) We define an L point sequence v n by the relation 

far*'" g n g M - 1 

Vn = \w-' L - n) ' n L-N + l^n<L. (18) 

arb itrary other n , if any 

If L is exactly equal to M + N — 1, the region in which v n is arbitrary 
will not exist. If the region does exist an obvious possibility is to in- 
crease M, the desired number of points of the z-transform we com- 
pute, until the region does not exist. 

Notice that v n could be cut into two with a cut between n = M — 1 
and n = L— N-\-V, if the two pieces were abutted differently, the result- 
ing sequence would be a slice out of the indefinite length sequence 
W~ n ' /2 . This is illustrated in Fig. 4. The sequence v n is defined the way 
it is in order to force the circular convolution to give us the desired 
numerical results of an ordinary convolution. 

(v) We compute the discrete Fourier transform of v„ and call it V r , 
r = Q,l,...,L-l. 

(vi) We multiply V r and Y r point by point, giving G> : 

G T = V T Y r , r = 0, 1, • • • , L - 1. 

(vii) We compute the L point inverse discrete Fourier transform 
0fc,of G>. 



Xn 




(a) 



(b) 



(c) 



yn 



N-l 



Y r 



N-l 



*W 



-n2/2 




L-i r 



(d) 



r-e- ARBITRARY- 




M-! 



L-N + i L-t n 



♦ Vr 




(f) 



o: 



(h) 



(L) 



L-l 




♦ 9k 



L-i r 



NOT USED 



M-1 



L-I k 



fX k 



M-1 



Fig. 4 — Schematic representation of the various sequences involved in the 
chirp z-transform algorithm: (a) input sequence x„ with N values, (b) weighted 
input sequence y n = A _ "TT'" 2 n x n ■ (c) discrete Fourier transform of y„. (d) values of 
the indefinite sequence W~ ni n . (e) sequence v n formed appropriately from segments 
of W~" 2 n . (f) discrete Fourier transform of v n . (g) product G> = Y r -V r . (h) inverse 
discrete" Fourier transform of G> . (i) desired M values of the z-transform. 



1258 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

(viii) We multiply g k by W k%/2 to give us the desired X k : 

X k = W k ' /2 g k , fc = 0, 1,2, ... ,M- 1. 

The ^ for fc ^ M are discarded. 

Figure 4 represents typical waveforms (magnitudes shown, phase 
omitted) involved in each step of the process. 

III. FINE POINTS OF THE COMPUTATION 

3.1 Operation Count and Timing Considerations 

An operation count can be made, roughly, from the eight steps just 
presented : 

(i) We assume that step i, that is, choosing L, is a negligible opera- 
tion. 

(ii) Forming y n from x n requires N complex multiplications, not 
counting the generation of the constants A~ n W n%n . The constants may 
be prestored, computed as needed, or generated recursively as needed. 
The recursive computation would require two complex multiplications 
per point. 

(Hi) An L point discrete Fourier transform requires a time /c F ft^ 
\0g2L for L a power of two, and a very simple fast Fourier transform 
program. More complicated (but faster) programs have more com- 
plicated computing time formulas. 

(iv), (v) The value of v n is computed for either M or N points, which- 
ever is greater. The symmetry in W~ n ' /2 permits the other values of v n 
to be obtained without computation. Again, v n can be computed re- 
cursively. The fast Fourier transform takes the same time as that in 
step Hi. If the same contour is used for many sets of data, V r need only 
be computed once, and stored. 

(vi) This step requires L complex multiplications. 

(vii) This is another fast Fourier transform and requires the same 
time as step Hi. 

(viii) This step requires M complex multiplications. 

As the number of samples of x n or X k grows large, the computation 
time for the chirp ^-transform grows asymptotically as something 
proportional to L log 2 L. This is the same sort of asymptotic depend- 
ence of the fast Fourier transform, but the constant of proportionality 
is bigger for the chirp 3-transform because two or three fast Fourier 
transforms are required instead of one, and because L is greater than 
N or M. Still, the chirp z-transform is faster than the direct com- 



CHIRP Z -TRANSFORM 1259 

putation of equation (10) even for relatively modest values of M 
and N of the order of 50. 

3.2 Reduction in Storage 

The chirp z-transform can be put into a more useful form for com- 
putation by redefining the substitution of equation (11) to read 

. (n - N ) 2 + k 2 - (fc - n + N.) 2 + 2N k 
nk- 2 

Equation (12) can now be rewritten as 

X„ = W h ' /2 W"' k 2! Xn A- n W tn - I, ' ) ' /a W- tk " ,+lr ' i ' / \ 

n-0 

The form of the new equation is similar to equation (12) in that the 
input data x n are preweighted by a complex sequence (A" n TF (n ~ Ar ° >,/2 ), 
convolved with a second sequence (W~ ln ~ N ' ) ' /2 ) i and postweighted by a 
third sequence (W k ' /2 W Nok ) to compute the output sequence X k . 
However, there are differences in the detailed procedures for realizing 
the chirp z-transform. The input data x n can be thought of as having 
been shifted by N samples to the left. For example, x„ is weighted by 
W N °' /2 instead of W°. The region over which W~ n ' must be formed, in 
order to obtain correct results from the convolution, is 

-N + I + N ^ n ^ M - 1 + N . 

By choosing N = (N — M)/2 it can be seen that the limits over which 
W~ n ' /2 is evaluated are symmetric; that is, W~ n ' /2 is a symmetric func- 
tion in both its real and imaginary parts. (Therefore, the transform of 
W~" 1/2 is also symmetric in both its real and imaginary parts.) It can 
be shown that by using this special value of N a only (L/2 + 1) points 
of W~ n ' /2 need be calculated and stored, and these (L/2 -f- 1) complex 
points can be transformed using an L/2 point transform*. Hence the 
total storage required for the transform of W~ "" is L + 2 locations. 
The only other modifications to the detailed procedures for evaluating 
the chirp z-transform presented in Section II are: 

(i) following the L point inverse discrete Fourier transform of step 

vii, the data of array g k must be rotated to the left by N locations, 

\ii) the weighting factor of the g k is W k ' /2 W Nck rather than W k ' /a . 

* The technique for transforming two symmetric L point sequences using one 
L/2 point fast fourier transform was demonstrated by J. Cooley at the fast 
Fourier transform workshop, Arden House, Harriman, New York, October 1968. 
The appendix summarizes this technique. 



1260 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

The additional factor W N ° k represents a data shift of N samples to 
the right, thus compensating for the initial shift and keeping the 
effective positions of the data invariant to the value of N used. 

Now we can estimate the storage required to perform the chirp z- 
transform. Assuming that the entire process is to take place in core, 
storage is required for V r which takes L + 2 locations, for y n , which 
takes 2L locations, and perhaps for some other quantities which we wish 
to save, such as the input or values of W*' /a or A~*W*' /a . 

3.3 Additional Considerations 

Since the chirp z-transform permits M ?* N, it is possible that occa- 
sions will arise where M ~S> N or N » M. In these cases, if the smaller 
number is small enough, the direct method of equation (10) is called 
for. However, if even the smaller number is large it may be appropriate 
to use the methods of sectioning described by Stockham. 3 Either the 
lap-save or lap-add methods may be used. Sectioning may also be used 
when problems, too big to be handled in core memory, arise. We have not 
actually encountered any of these problems and have not programmed 
the chirp z-transform with provision for sectioning. 

Since the contour for the chirp z-transform is a straight line seg- 
ment in the s-plane, it is apparent that repeated application of the 
chirp z-transform can compute the z-transform along a contour which 
is piecewise spiral in the z-plane or piecewise linear in the s-plane. 

Let us briefly consider the chirp z-transform algorithm for the case 
of z k all on the unit circle. This means that the z-transform is like a 
Fourier transform. Unlike the discrete Fourier transform, which by 
definition gives N points of transform for N points of data, the chirp 
z-transform does not require M = N. Furthermore the z k need not 
stretch over the entire unit circle but can be equally spaced along an 
arc. Let us assume, however, that we are really interested in comput- 
ing the N point discrete Fourier transform of N data points. Still the 
chirp z-transform permits us to choose any value of N, highly com- 
posite, somewhat composite, or even prime, without strongly affecting 
the computation time. An important application of the chirp z-trans- 
form may be computing discrete Fourier transforms when N is not a 
power of two and when the program or special purpose device avail- 
able for computing discrete Fourier transforms by fast Fourier trans- 
form is limited to when N is a power of two. 

There is also no reason why the chirp z-transform cannot be ex- 
tended to the case of transforms in two or more dimensions with simi- 
lar considerations. The two dimensional discrete Fourier transform 



CHIRP Z-TRANSFORM 1261 

becomes a two dimensional convolution which can be computed by- 
fast Fourier transform techniques. 

Caution: For the ordinary fast Fourier transform the starting point 
of the contour is still arbitrary; merely multiply the waveform x n by 
A~ n before using the fast Fourier transform and the first point on the 
contour is effectively moved from z = 1 to z = A. However, the con- 
tour is still restricted to a circle concentric with the origin. The angular 
spacing of z k for the fast Fourier transform can also be controlled to 
some extent by appending zeros to the end of x n before computing the 
discrete Fourier transform (to decrease the angular spacing of the 
z k ) or by choosing only P of the N points x n and adding together all 
the x n for which the n are congruent modulo P; that is, wrapping the 
waveform around a cylinder and adding together the pieces which 
overlap (to increase the angular spacing) . 

IV. APPLICATIONS OF THE ALGORITHM 

Because of its flexibility, the chirp 2-transform algorithm discussed 
in the Section III has many potential applications. 

4.1 Enhancement of Poles 

One advantage of the chirp 2-transform algorithm over the fast 
Fourier transform is its ability to evaluate the 2-transform at points 
both inside and outside the unit circle. This is important in the investi- 
gation of systems whose transfer functions can be represented as ratios 
of polynomials in z; that is, in finding poles and zeros of a linear sys- 
tem. By evaluating the transform off the unit circle, one can make the 
contour pass closer to the poles and zeros of the system, thus effectively 
reducing the bandwidths and sharpening the transfer function. 

For example, a five-pole system was simulated at a 10 kHz sam- 
pling frequency. The poles were located at center frequencies of 270, 
2290, 3010, 3500 and 4500 Hz with bandwidths of 30, 50, 60, 87 and 
140 Hz, respectively. The z-plane pole positions are shown in Fig. 5. 
(Those familiar with speech will recognize these numbers as resonance 
positions appropriate for the vowel i in the word beet.) The input to 
the system was a periodic impulse train of period 100 samples; that is, 
100 pulses per second. Impulse invariant techniques were used to 
simulate the system. 6 The 2-transform of one period of steady state 
data (100 samples) was evaluated on two spirals outside the unit 
circle, one on the unit circle, and two spirals inside the unit circle. 
Figure 6 shows the five contours as they would appear in the s-plane 
and the s-plane pole positions. The contours are seen to be equi- 



1262 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

Z- PLANE 




Fig. 5 — Representation of the z-plane locations of the poles of the linear sys- 
tem simulated in the text. 

angularly spaced. The five sets of magnitude curves are shown in Fig. 
7. The transform was evaluated at 50 equally spaced points from 
to 4900 Hz, corresponding to <p = — 1/100. The sharpening of the 
magnitude response in the region of the poles is quite pronounced. 
Figure 6 indicates that contour 5 should be near optimum since it 
intersects three of the poles. 

This example is a somewhat idealized case in that spectral samples 
were taken every 100 Hz; that is, at the harmonics of the funda- 
mental frequency. Figure 8 shows the case for spectral data taken 
every 25 Hz on contour 5 of Fig. 6, along with the case where the 
spectral resolution is the same as shown in Fig. 7. This figure places 
in evidence the true nature of the z-transform of a finite number of 
samples. It is clear from equation (2) that X(z) has no poles any- 
where in the z-plane except at z = 0. There are instead N—l zeros 
which manifest themselves in the ripples seen in the upper curve of 
Fig. 8. In many cases the poles of the original system which generated 
the samples are still in evidence because the zeros tend to be arrayed 
at approximately equal angular increments except at the locations of 
the original poles. Hence a pole usually manifests itself by an absence 
of zeros in the vicinity of that pole in the z-plane. Zeros of transmis- 
sion are often masked by these effects when only a finite number of 
samples are transformed. Examples of this effect are given after equa- 
tion (23). 

It is of interest to examine the ability of the chirp z-transform 
algorithm to determine the bandwidth of a pole as well as its center 



CHIRP Z -TRANSFORM 



1263 



frequency. To investigate this point, synthetic samples were generated 
with the bandwidth of the lowest pole (Bx) variable from 10 to 320 
Hz by factors of 2; all other bandwidths and center frequencies were 
held at values used in the previous example. Again a fixed 100 pulse 
per second source excited each of the systems. Figure 9 illustrates the 
six sets of poles and three contours used in the investigation. Contour 
3 extends into the right half-plane (spiral outside the unit circle) and 
is only close to the lowest pole. Contour 2 corresponds to the unit circle 
in the z-plane (that is, the discrete Fourier transform). Contour 1 is 
an appropriate left-half plane contour (spiral inside the unit circle) 
used for investigating this system. The resulting set of 18 magnitude 
curves (six different sets of poles and three contours) are shown in 
Fig. 10. The rows of Fig. 10 show magnitude curves with a fixed band- 
width and variable contour, whereas the columns show curves on the 
same contour with variable bandwidths. There are 801 points plotted 
in each curve in the range to 5,000 Hz. 

Looking down any column it is seen that as B± increases, the level 
of the first resonance decreases steadily. The variation in fine spectral 
detail resulting from the distribution of zeros in the neighborhood of 
the poles of the original system is seen clearly in column 1. For ex- 



4000 



X 

z 3000 \- 



O 2000 

s 



1000 



LEFT HALF-PLANE 




-150 



-50 50 100 

FREQUENCY IN HZ 



Fig. 6 — The s-plane locations of the poles of Fig. 5 and five contours for 
evaluation of the z-transform. 



1264 



THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1909 




1000 2000 3000 4000 5000 
FREQUENCY IN HZ 



Fig. 7 — Magnitude curves corresponding to evaluation of the z-transform on 
the five contours of Fig. 6. 




1000 2000 3000 4000 5000 

FREQUENCY IN HZ 



Fig. 8 — A comparison of high resolution and low resolution evaluations of 
the z-transform. The spacing of points is 25 Hz in the upper curve and 100 Hz 
in the lower curve. 



CHIRP Z-TRANSFORM 



1265 



4000 



a 2000 - 




-320 



-50 

FREQUENCY IN HZ 



150 



Fig. 9 — The pole locations and contours used to investigate the possibility of 
bandwidth determination using the chirp z-transform. 



ample, the fifth resonance at 4,500 Hz is difficult to find in the upper 
plots and almost missing in the lower plots, because of the presence of 
a zero at the pole position. Furthermore, the frequency at which the 
magnitude is minimum, that is, the closest zero to contour 1 shifts 
from 2,500 to 2,700 to 800 to 1,100 Hz as bandwidth increases. 

The plots in columns 2 and 3 show little or no variation from about 
2,000 to 5,000 Hz where the appropriate contours are generally far 
away from the zeros of the distributions. The resonance at 4,500 Hz 
is always easy to locate on these plots, thus indicating the desirability 
of both detailed close-up examination of the transform (as on con- 
tour 1) and less detailed, further away looks at the magnitude curve 
(as on contours 2 and 3). The magnitude curves in the regions to 
2,000 Hz are still fairly sensitive to the exact zero distribution for 
contour 2, and slightly sensitive for contour 3. It would appear from 
Fig. 10 that there are cases when bandwidth can be determined either 
from the width or the magnitude of the resonance. Further investiga- 
tion is necessary before quantitative techniques for determining band- 
widths are available. 

The choice of the optimum contours is highly dependent on the 
locations of the poles of the original system. In general there is no 



1266 



THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 



I 


N 

X 




l\l 


II 


II 


DO 


CO 




snaanaa ni aaruiNsvw 



CHIRP Z -TRANSFORM 



1267 



single contour on which all the poles are located since these contours 
are essentially lines of constant Q = (center frequency) /(bandwidth). 
Hence the choice of contour is highly dependent on which of the sys- 
tem poles is of most interest in the particular problem. However, some 
interesting observations can be made from studying magnitude curves 
for systems whose poles are constant Q poles. Such a system was 
simulated by keeping the pole center frequencies at the values used 
previously and setting the Q of each of the poles to 20. A 100 pulse 
per second impulse train was again used to excite the system and one 
period of steady state data was analyzed along contours 1 to 7 shown 
in Fig. 11. The pole positions of the system are shown in this figure 
and are seen to coincide with contour 6 exactly. The magnitude curves 
for these seven contours are shown in Fig. 12. These are high resolution 
spectra containing 801 points from to 5,000 Hz. Notice that both 
magnitude curves 5 and 6 accentuate the poles equally except in the 
region of the fifth pole where curve 5 appears slightly better than 
curve 6. The fact that this occurs is not surprising in view of the fact 
that the pole is really manifested by an absence of a zero in an array 
of zeros of approximately the same magnitude and angular spacing. 
Another anomaly which can be attributed to the way in which the 



4000 - 



2 2000 - 




-200 -150 -100 

FREQUENCY IN HZ 



Fig. 11 — The s-plane locations of constant Q poles and contours on which the 
z-transform was evaluated. 



1268 



THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 




90 -, 
70 

50 
30 h 




105 

85 
65 

45 





1000 2000 3000 4000 
FREQUENCY IN HZ 



Fig. 12 — Magnitude curves for contours of Fig. 11 (constant Q poles). 

zeros are distributed is evidenced by comparing curves 5 and 7. Based 
on the relative positions of the two contours with respect to the poles 
we would expect the magnitude curves for these contours to be identi- 
cal, but the comparison shows that this is not the case in actual com- 
putation. This is the result of the fact that the zero distribution is not 
exactly symmetric so that contours which pass very close to the zeros 
may look considerably different from one another. 

A final point of interest in Fig. 12 is the linear component in the 
last three curves which dominates at high frequencies. This effect is 
also shown in Fig. 14. Figure 13 shows the five contours used in ob- 



CHIRP Z -TRANSFORM 



1209 



taining the log magnitude plots in Fig. 14. It is clear that when the 
contour passes inside the original pole locations (and therefore inside 
the array of zeros in the 2-plane) , the log magnitude function exhibits 
a definite linear component. This effect is easily explained when X(z) 
is written in the form 



X(z) - Dz-' N ~ l) 11(1 -a7 l z). 



(19) 



where the a r 's are the zeros of X(z). If we evaluate equation (19) at 
z = W~ k exp (—j2TT(p k) we obtain for the magnitude 

\X k \= \D\ WV' l} II 111 - a; l W: h exp i-frrJk)]\- (20) 
For plotting in dB we define 
20 log,„ | X k | = 20 log.o \D\ + 20(N - l)k log 10 W 

+ Z 20 log 10 |[1 - a: x W: k exp (-j2ir^.fc)]|. (21) 

r-l 

In the examples we have shown, almost all of the zeros a r have mag- 
nitudes slightly less than 1. Thus for contours inside these zeros (cor- 



9, 2000 




-600 -500 -400 -300 -200 

FREQUENCY IN HZ 



Fig. 13 — Contours and pole locations used to study the effect of passing in- 
side the pole locations. 



1270 



THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 



85 


A A A . 


70 
55 


A MV/V\ 




i i i i 










110 




90 




70 






^ — r^ i i i 










140 




no 




80 




60 


■ / ^— ■ -r i i i 


240 


^^^ 


180 


^ — ^^ 


120 


^^~^ 


60 


— -<^ i i i 


280 




200 




120 




40 


•* i i i i 



NO. 2 



NO. 3 



NO. 4 



1000 2000 3000 4000 5000 
FREQUENCY IN HZ 

Fig. 14 — Magnitude curves for the contours of Fig. 13 showing a large linear 
component resulting from the N — 1 poles of X(z) at z = 0. 

responding to contours 2 through 5 in Fig. 13), W is greater than 1. 
Thus each term in the sum on the right side of equation (21) tends to 
decrease as k gets larger. In contrast, the second term on the right 
side of equation (21) represents a linear component with slope equal 
to20 (iV— 1) log 10 F . 

In Fig. 14, the rath curve corresponds to a value of W = e 47rm/10 '° 00 , 
m = 1, 2, 3, 4, 5. The value of <p is - 1/100, N is 100, and the sam- 
pling rate is 10 kHz. Thus the frequency going from to 5 kHz cor- 
responds to k going from to 50. For example, in the fifth curve W = 
e 2 0w /io.ooo anc j t h e s i p e should be 

20(tf - 1) log.. W. = ^y" - 5.4. (22) 

Thus the total dB change in going from to 5 kHz should be on the 
order of 50(5.4) = 270 dB. In Fig. 15a we show this case again. In 
Fig. 15b we show the result of evaluating 



2; w.r = 2 wr " I, Ki 



(23) 



»— (AT-l) 



CHIRP Z-TRANSFORM 



1271 



using the same value of W and <p . Notice that this should remove the 
second term in equation (21) leaving the other terms unaffected. This 
observation is substantiated by Fig. 15b since the value at 5 kHz is 
very nearly 270 dB less than in Fig. 15a. Notice also that some of the 
resonances are still in evidence although not as clearly defined because 
the contour is passed relatively far inside the zeros of X(z). 

Another interesting question was investigated using this technique. 
As shown above, when the response of a linear system is truncated by 
repetitively pulsing the system and transforming a finite number of 
samples, the z-transform has only zeros except for poles at z = 0. 
However, the poles of the original system function can still be located 
in magnitude response curves by the absence of zeros in the ap- 
propriate regions. The question arises about what happens when the 
system function contains zeros. Suppose h (nT) is the impulse response 
of a linear system with both poles and zeros in its z-transform H(z). 
Since H(z) has poles, h(nT) will be an infinite sequence. There is 
clearly no reason to expect that the transform of only N of these 
samples will have zeros at the same location as the zeros of H(z). 
However, the system zeros can be expected to have an effect on the 
distribution of zeros of the truncated z-transform. 

To illustrate this point, the system response used in the above ex- 
amples was modified by passing the output waveform through a sys- 
tem whose transfer function consisted of a complex conjugate zero 
pair. A periodic 100 pulse per second source was again used to excite 
the system and one period of steady state data was analyzed. The 
system pole-zero pattern and the contours of analysis are shown in 




I0OO 2000 3000 4000 5000 
FREQUENCY IN HZ 



Fig. 15 — Magnitude curves obtained by evaluation of the z-transform on 
contour 5 of Fig. 13: (a) with the effect of the N — 1 poles at z = 0; (b) with 
the N — 1 poles removed by shifting the sequence x» by N positions to the 
left. 



1272 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

Fig. 16. In one simulation a zero was placed at point A (500 Hz, 12.5 
Hz), and in a second simulation the zero was at point B, (2,500 Hz, 
60 Hz). The analysis was made at 801 points from to 5,000 Hz along 
these contours. The resulting magnitude curves, along with a set of low 
resolution curves where the magnitude was computed every 100 Hz 
from to 5,000 Hz, are shown in Figs. 17 and 18. The data of Fig. 17 
are for the case where the transmission zero was at 500 Hz whereas 
for Fig. 18 the zero was at 2,500 Hz. 

The high resolution data of Fig. 17 show no strong indication of 
the transmission zero; whereas the transmission poles are still very 
much in evidence. The low resolution data (evaluated at harmonics 
of the source) does indicate the presence of a zero along contour 1, 
but along the other contours the case is not so clear. The most unusual 
observation is that along contour 3, the contour closest to both the 
transmission zero and the poles, there is little or no indication of the 
zero; whereas the poles are still strongly in evidence. Along contour 
4, at the high frequency end, there is noise in the magnitude spectrum. 
The source of this noise is discussed in Section V. 

The indications from Fig. 17 are that a transmission zero can be 
more easily located on contours which are far from the zero than on 



5000 



4000 - 



9, 2000 - 



1000 - 




-120 -90 -60 
FREQUENCY IN HZ 



Fig. 16 — The s-plane locations of poles and zeros (at A and B) and contours 
used in studying the effect of zeroes on the magnitude curves. 



CHIRP Z-TRANSFORM 



1273 







HIGH RESOLUTION DATA 
(801 POINTS) 


190 






160 






130 






100 




" 1 1 1 1 








190 




A 7 Nw/ V_w^w 


160 






130 


^ 


fvvttftnr 






i • i i i 








190 
160 




w-v- 


130 


*y 


««*^ 






1 1 1 1 


280 










240 






200 






160 
IPO 




-~"l 1 1 ! 



LOW RESOLUTION DATA 
(51 POINTS) 



1000 2000 3000 4000 5000 



FREQUENCY IN HZ 




NO.? 



N0.4 



5000 



Fig. 17 — Magnitude curves for a zero at 500 Hz (position A in Fig. 16). 

contours which traverse it. Furthermore it is much easier to locate on 
a low resolution spectrum than on a high resolution spectrum. Hence 
zeros, unlike poles, are not generally easy to locate from spectra. 

The zero of Fig. 17 was at 500 Hz and in a region where the high 
resolution spectra displayed a large amount of ripple from the trunca- 
tion zeros of the data. Figure 18 shows similar magnitude curves for 
the zero at 2,500 Hz, a region with much less ripple in the spectrum. 
The magnitude response curves show effects entirely similar to those 
of Fig. 17. The zero is most easily beatable for contour 1, the stand- 
ard fast Fourier transform. In contour 3, which again passes through 
the zero, there is no indication of the zero. Also the low resolution data 
tends to show the zero better than the high resolution data. One im- 
portant implication of these results is that one could not use these 
techniques to accurately find the position of complex transmission 
zeros. In many cases it would be difficult to differentiate between dips 
in the spectrum between poles and dips caused by complex zeros, thus 
indicating the difficulty of locating even the center frequency of a 
zero. 

The chirp z-transform algorithm has been applied to the spectral 
analysis of speech in order to aid in automatic detection of the time 
varying resonances (poles or formants) of speech. Voiced speech can 



1274 



THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 



HIGH RESOLUTION DATA 
(801 POINTS) 



LOW RESOLUTION DATA 
(51 POINTS) 




N0.1 



N0.2 



NO.3 



N0.4 



1000 2000 3000 4000 5000 1000 

FREQUENCY IN HZ 



2000 3000 4000 5000 



Fig. 18 — Magnitude curves for a zero at 2,500 Hz (position B in Fig. 16). 

be modelled as the convolution of a source waveform with a vocal 
tract impulse response. The vocal tract impulse response is essentially 
a sum of damped exponentials, each exponential corresponding to a 
mode or pole of the vocal tract transfer function. It is of interest to 
speech researchers to detect these time varying resonances. The chirp 
2-transform algorithm has been applied to individual periods of voiced 
speech with a high degree of success. Figure 19 shows the result of ap- 
plying the chirp 2-transform algorithm along the two contours shown 
at the upper left of the figure, to a period of voiced speech. The upper 
contour corresponds to the standard fast Fourier transform contour; 
the lower to a suitably chosen spiral contour. The magnitude function 
along the upper contour indicates a single wide peak in the region 
2,000 to 2,500 Hz, whereas the magnitude along the lower contour 
shows two isolated peaks in this region corresponding to the physical 
knowledge that there actually are supposed to be two peaks in this 
region. Variations on the chirp 2-transform algorithm for spectral 
analyses of speech have been studied and will be reported on in a 
subsequent paper. 7 

4.2 High Resolution, Narrow Band Frequency Analysis 

One very useful application of the chirp z-transform algorithm is the 
ability to efficiently evaluate high resolution, narrow frequency band 



CHIRP Z-TRANSFORM 



1275 



spectra. Using standard fast Fourier transform techniques, in order to 
achieve a frequency resolution of ^ AF, with a sampling frequency of 
the data of l/T, requires N ^ 1/(T-AF) points. For very small AF, 
this implies very large values of N. The crucial issue is that what is 
often required is high resolution for a limited range of frequencies and 
low resolution for the remainder of the spectrum. An example of such a 
circumstance is the design of band-pass or low-pass filters. Usually 
what is desired is a microscopic look at details of the frequency response 
in the pass-band and only a gross look outside the pass-band. 

The chirp 2-transform algorithm is extremely well suited for such 
cases since it allows selection of initial frequency and frequency spac- 
ing, independent of the number of time samples. Hence high resolu- 
tion data over a narrow frequency range can be attained at low cost. 

To illustrate these points, simple rectangular band-pass niters were 
simulated by symmetrically truncating a delayed impulse response. 
The impulse response used was 



h(nT) = a sin [ir(F 2 - F,)(n -\- m)T] 
•cos [t(F 2 + F^n - \ - m)T), 



< n < 2ra 



(24) 



S-PLANE 





3500 4000 



Fig. 19 — Magnitude curves from evaluation of the z-transform of one period 
of natural speech. The contour for the upper plot is the unit circle in the 2-plane 
while the contour for the lower curve is a spiral inside the unit circle. 



1276 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

where 

2ra = number of terms in the truncated impulse response 
1/T = sampling frequency = 10,000 Hz 
F x = lower cutoff frequency in Hz 
F 2 = upper cutoff frequency in Hz. 

Values for m of 100 and 500 were used with F x = 900 Hz and F 2 = 
1,100 Hz. Figure 20 shows plots of equation (24) for these two cases. A 
standard 1,600 point fast Fourier transform was calculated and the 
magnitude response for m = 100 is shown in the upper half of Fig. 21. 
In order to investigate the pass-band and transition region more care- 
fully the chirp z-transform algorithm was used to give a 1.25 Hz 
resolution over the band from 500 to 1,500 Hz. The contour used was 
identical to the contour for the fast Fourier transform. The resulting 
magnitude response curve is shown in the lower half of Fig. 21. To 
achieve this high a resolution would have required an 8,000 point 
fast Fourier transform, instead of the 1,000 point transforms actually 
used. (Similar expansions of regions of the phase curve were made for 
this filter but are not shown.) 

Figure 22 shows similar effects for the case m = 500. The applica- 
bility of the chirp z-transform algorithm for such frequency expan- 
sions is a powerful tool for close examination of small frequency 
bands, as well as for debugging implementations of digital filters. 
For example, one could easily check if a desired filter met its design 
specification of in-band ripple, transition ratio, and so on. 8 

One situation where the chirp z-transform algorithm may be quite 
useful is when we are confronted with an extremely long sequence for 
which we desire a fine grained spectrum over a narrow band of fre- 
quencies. Suppose we have a sequence of P samples and desire M 
spectral samples where M <<C P. That is, we wish to evaluate 

p-\ 
X k = £ x n A~ n W k , k = 0, 1, ••• ,M - 1. (25) 

n-0 

The sum in equation (25) can be broken up into r sums over N points 
as follows 



X k m 2 A-" V W 



kqX 



2 x n+QlV A- n W nk , k - 0, 1, ■ • • , M - 1 

n-0 J 

(26) 



where rN ^ P. Each of the r sums in the brackets can be evaluated 
using the chirp z-transform algorithm, requiring storage on the order of 



CHIRP Z-TRANSFORM 



1277 



10000 




4^4^* 



300 400 500 600 
SAMPLE NUMBER 



700 800 900 1000 



Fig. 20 — The impulse responses of simple band-pass filters. 



3(N + M — 1) locations. In addition we require 2M locations in which 
to accumulate the M complex values of the transform. Although 2 
fast Fourier transforms and 2M complex multiplications are required 
for each of the r transforms, it is quite possible that a saving in total 
time may result from this method as opposed to evaluation of a P 
point transform using auxiliary storage such as drum, disk or tape. 



1278 



THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 



160 


- A 


100 


120 


Mgf^ ^mmrnmm- 




80 


\\\\W m]\\mm 


pPTOfm 


PfRfftam^ 


in 

i 4 ° 

o 

Q 


i i i i i 




i i 


i 



2000 3000 

FREQUENCY IN HZ 




900 1000 1100 

FREQUENCY IN HZ 



Fig. 21 — Frequency response curves for upper impulse response (200 samples) 
in Fig. 20. Upper curve obtained with 1,600 point fast Fourier transform (resolu- 
tion 6.25 Hz). Lower curve obtained with chirp z-transform algorithm (1.25 Hz 
resolution) . 



4.3 Time Interpolation or Sampling Rate Changing 

The flexibility of the chirp ^-transform algorithm for obtaining high 
resolution in frequency has been explained and illustrated in Section 
4.2. A similar procedure applies to interpolation between samples of a 
bandlimited time function using samples of the frequency spectrum. 9 
In this section, we discuss how the discrete Fourier transform can be 
used to perform interpolation on a set of samples and the advantages 
and disadvantages of using the chirp 2-transform algorithm for this. 

4.3.1 Bandlimited Interpolation Using the discrete Fourier transform 
Assume that we have available N samples x(nT) , n = 0, 1, 2, . . . , 
N — 1, of a bandlimited waveform x{t). The sampling interval T is 
assumed less than or equal to the Nyquist interval. The total time in- 
terval spanned by these samples is therefore NT seconds. We wish to 
obtain equally spaced samples of x(t) at a sampling interval T f , where 



CHIRP Z -TRANSFORM 



1279 



T is less than or equal to the Nyquist interval. These samples are 
denoted by x(mT') , m = 0, 1, . . . , N' - 1, where N'T' = NT. (Notice 
that we are assuming N' is an integer. This assumption will be dropped 
later.) 

If all the samples x(nT) are available, the samples x{mT') can be 
obtained from 



x(mT') = £ x ( nT ) 



sin ^ (mT r - nT) 
| (mT' - nT) 



(27) 



Thus the interpolation can be viewed as the result of convolving the 
interpolation function [sin (7rt/T)]/[(ir/T) t] with the samples re (nT) 
and then resampling with period T'. It is well known that convolution 
may be done using the discrete Fourier transform, and we will show 
how the resampling can also be affected by properly augmenting the 



160 


- n 


m= soo 


120 


mmmdtW \\ 




80 


Wnvfv "\ 


TWnpfyw^^ 


m 

w 40 

u 

* n 


i i 


M ' II ' l ' | j 1 ' 

i i i i i i i 



2000 3000 

FREQUENCY IN HZ 



5000 




900 1000 1100 

FREQUENCY IN HZ 



Fig. 22 — Frequency response curves for lower impulse response (1,000 sam- 
ples) of Fig. 20. Upper curve obtained with 1,600 point fast Fourier transform 
(resolution 6.25 Hz). Lower curve obtained with chirp z-transform algorithm 
(1.25 Hz resolution). 



1280 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

transform with zeros. Because the discrete Fourier transform uses only 
a finite number of samples, we shall encounter errors similar (but not 
identical) to using only N terms in equation (27) . 

The discrete Fourier transform of the given samples is 



XM = g x(nT) exp (- j || knTJ 



k = 0, 1, ••• ,N - 1. (28) 

(Notice that we have changed notation in this section in order to make 
explicit the number of samples and the sampling period.) We define 

X&k) = X N {k)H N (k), (29) 

where H N (k) is the N point discrete Fourier transform of the inter- 
polation function to be convolved with the samples x(nT). [Notice 
that this convolution is equivalent to cyclic convolution of a periodic 
impulse response h(nT) with the samples x{nT) .] 

In order to change the sampling to period T', we split X/,{k) about 
k = N/2 and expand (by inserting zeros) or contract (by discarding 
zeros) the transform according to the following equations 



X' N .(k) = xm 


10 ^ k < N'/2 
lo ^ k < N/2 




N' <N 

N' > N 


(30a) 


= 


k = N'/2 




N' < N 


(30b) 


= 2Aat(/c) 


k = N/2 




N' > N 


(30c) 


= 


N/2 <k <N' - N/2 




N' > N 


(30d) 


= \X' N (k - 


- N' + N) k = N' - 


N/2 


N' < N 


(30e) 


= X' N (k - 


N , + J N>/2 <k<N> 
IN' - N/2 < k < N' 


N' <N^ 
N' > N 


(30f) 



Equations (30b) , (30c) , and (30e) are required only when N' and N 
are even integers and equation (30d) is required only when N' > N. 

The N' point inverse discrete Fourier transform of X' N , (k) is defined 
to be 

x'{mT) = ~ g XUk) exp (; ^ kmT') , 

m = 0, 1, ••• ,N f - 1. (31) 



CHIRP Z-TRANSFORM 1281 

For example, if N is even and N' > N, we can show using equations 
(28), (29), and (30) that 

I AT/2 Tat-i / 9 ^ 

x'(mT') = p k S *»<*)[ § *(») ex P (-J ^f ^ 

■exp(;^7/fcmr), (32) 

where m = 0, 1, • • • , N' — 1, and the terms corresponding to k = ±N/2 
are understood to be multiplied by 1/2 since N is even. By interchanging 
the order of summation and using the fact that N'T' = NT, we obtain 

x'imT') = tt/ E x(nT)h(mT' - nT) (33) 

*V „-0 

where 

h(mT' ~nT) = jj ^ H„(k) exp j f^ k(mT' - nT) J ■ (34) 

Notice that equation (33) has the desired form for interpolation [see 
equation (27)], however the values of x'imT') are clearly not exactly 
equal to the desired interpolated values x{mT / ), This is so because 
only N samples are used and because of the form of h (mT — nT) . As 
an example, suppose 

HM = fl -N/2 < k < 2V/2. (35) 

li k= ±N/2 

(This is equivalent to splitting X N {k) at k = N/2 and inserting N' — 
N zeros between the two halves of the transform) . If we evaluate 
equation (34) for this case, we obtain 

sin ~ (mT' - nT) 

h(mT' - nT) = (36) 

Ntanj^(mT'-nT) 

This function is plotted in Fig. 23a where = (mT' — nT)/T, and 
N = 8. Clearly 

sin ~ (mT' - nT) 

h(mT' - nT) « (37) 

| (m2" - nT) 



1282 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

(a) 



N = 8 




K 





t ^^°~~- 




(b) 










1 I 


1 1 


1 


1 




T 


2T 


1 1 




NT 

I 


t 


1 


2 






N 


n 




N' m 



Fig. 23 — An illustration of bandlimited interpolation using the discrete Fourier 
transform: (a) the periodic function which is convolved with the original sam- 
ples; (b) a bandlimited time function showing samples with spacing T; (c) How 
the interpolated values xfimT') are formed. 

when 7r(mT' — nT)/NT\s small, that is, in a region where h(mT' — nT) 
is significantly different from zero. Figure 23b shows a segment of a 
waveform x{t) and samples x(nT). In Fig. 23c, we have shown just 
two of the terms in equation (33). This figure places in evidence the 
nature of the interpolation which is performed. The errors are likely 
to be greatest at either end of the segment, since the interpolated values 
at one end depend on the samples at the other end in a way which is 
not at all consistent with equation (27). The error caused by this 
effect will be most significant in the regions ^ mT' < 2T and 



CHIRP Z-TRANSFORM 



1283 



T(N — 2) ^ mT' ^ TN. The remainder of the values will have es- 
sentially the same error associated with using only N terms in equa- 
tion (27). 

Notice that equation (36) is not the only interpolation function 
which can be used. Other choices of H N (k) may lead to interpolation 
functions which are in some sense more desirable. For example, Fig. 
24 shows four different choices for H N (k) and their associated inter- 
polation functions or impulse responses. (The impulse responses were 
shifted modulo N' as an aid in plotting.) It can be seen from Fig. 24, 
that removing the sharp cutoff in H N {k) greatly shortens the effective 
duration of the impulse response, thus tending to minimize the end 
effects discussed previously. Clearly the approximation to (sin wt)/vt 
interpolation is not as good as equation (36) , but in many cases such 
smoothing of the interpolated values may not be objectionable. 



1.0 
0.8 

o.e 

0.4 

0.2 




800 

w 600 

3 400 

£ 200 



-200 



_ IMPULSE 
RESPONSE 



FREQUENCY 
RESPONSE 



1000 2000 3000 4000 
FREQUENCY IN Hz 



20 30 40 

TIME IN SAMPLES 

(a) 



1.0 

0.8 
0.6 
0.4 
0.2 



800 

u 600 
400 

( 200 


-200 



FREQUENCY 
RESPONSE 



1000 2000 3000 4000 5000 
FREQUENCY IN HZ 



IMPULSE 
RESPONSE 



40 60 80 100 120 140 160 
TIME IN SAMPLES 
(b) 



1.0 
0.8 
0.6 
0.4 

0.2 




800 

600 
UJ 

3 400 

£ 200 


-200 



FREQUENCY 
RESPONSE 



1000 2000 3000 4000 
FREQUENCY IN HZ 



IMPULSE 
RESPONSE 




20 30 40 

TIME IN SAMPLES 

(O 



i.o 

0.8 
0.6 
0.4 

0.2 




800 

600 
J 

3 400 

C 200 



-200 



FREQUENCY 
RESPONSE 



1000 2000 3000 4000 5000 
FREQUENCY IN HZ 



IMPULSE 
RESPONSE 



20 40 60 80 100 120 140 160 
TIME IN SAMPLES 

(d) 



Fig. 24 — A set of four simple frequency responses and corresponding im- 
pulse responses which could be used for interpolation. 



1284 



THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 



4.3.2 Computational Considerations in Bandlimited Interpolation 
In Section 4.3.1 we discussed a method of bandlimited interpolation 
based on the discrete Fourier transform. The operations involved are 
summarized in Fig. 25. The sequence {X N {k) } may be evaluated using 
the fast Fourier transform. In this case, N will be restricted to a value 
compatible with the available fast Fourier transform routine, for ex- 
ample, N would be highly composite. The transform is then multipled 
by H N (k) and expanded or contracted according to equation (30). 
Then we must compute the inverse discrete Fourier transform with N' 
points. This can be done using the fast Fourier transform provided 
that 

(i) JV' = NT/1" is an integer and compatible with the available fast 
Fourier transform routine. 

(ii) Enough high speed storage (a minimum of N' locations) is 
available. 

[Notice that i applies for either N' > N or N' < N while ii will probably 
not be a problem except when N' y> N.] 

In many cases it may not be possible to meet one or both of the 
above conditions; then the chirp z-transform algorithm can be very 
useful. The N' point inverse discrete Fourier transform may be com- 
puted using W = 1 and <p = + l/iV', where N' need not even be an 
integer. Thus we can compute M interpolated values using 



N/2 

2 

k—N/2 



1 \LLz 

x'(mT) = ^7 W m ' n E [X' N .(k)W" /2 ]W- ( 



-(m-*)V2 



where X£.(k) is determined by equation (30) and 
X' K .(-k) = X' N .(N' ~k), k = l,2, 



,N/2. 



(38) 



(39) 



x(nT) 



DISCRETE 

FOURIER 

TRANSFORM 



H N (k) 



(Xhr 



x N (k) v-^ X ' N (k) 



n = o,i,...,N-i 



N' 



EXPAND 

OR 

CONTRACT 



xW(k) 



N' 



INVERSE 

DISCRETE 

FOURIER 

TRANSFORM 



x'(mT') 



m = o,i,...,N-i 



Fig. 25 — Illustration of the steps involved in bandlimited interpolation using 
the discrete Fourier transform. 



CHIRP Z-TRANSFORM 1285 

Assuming that the transform of W~ k ' /2 is available, equation (38) 
can be evaluated using two L point fast Fourier transforms where L 
is the smallest integer which is greater than N + M — 1 and which is 
compatible with an available fast Fourier transform routine. 
Alternatively we can evaluate 

»W = W *■"■ 2! [Y'Ak)W" /2 ]W- l - k) ' /2 (40) 

1 v ,. = 

where 

Y> N .(k) = \ X ' Ak) k = ° • (41) 

[2Xit.(k) < k ^ N/2 

It can be shown that 

y'{mT') = x'{mT') - jx'(mT') (42) 

where x'(mT') is the inverse discrete Fourier transform of 

l' N .{k) = jsgn„.(k)-X' N .(k) (43) 

and 

k = 0, N'/2 

sgn jV - (fc) - I 1 < k < N'/2. (44) 

-1 N'/2 < k < N' - 1 

From equation (41) and (44) it can be shown that x~'(mT') is an approxi- 
mation to the Hilbert transform of x'{mT). In this case we require at 
least {N/2 + M — 1) point transforms to compute M interpolated 
values. This is at the expense of not being able to do two interpolations 
at once as is possible with equation (38) (obtaining one interpolation 
as a real output and one as an imaginary output); however we do obtain 
an approximation to the Hilbert transform of x'iinT') which may be of 
value in some applications. 

If sufficient core storage is not available to compute an N' point 
fast Fourier transform, we can compute the interpolated values in 
sections and piece these sections together as is commonly done in high 
speed convolution. The chirp ^-transform algorithm allows us to com- 
pute as many as 2M interpolated points at a time, where M can be 
chosen so that the fast Fourier transforms can be done using only core 
storage. Probably the most significant advantage, though, is the 
ability to efficiently interpolate to arbitrary sampling intervals. 



128G 



THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 



As an example of the ideas discussed in this section, consider the 
waveforms shown in Fig. 26. Figure 26a shows 500 samples of a speech 
waveform where the sampling rate was 20 kHz. (T = 5 X 10~ B sec- 
ond). The samples are connected by straight lines in the figure. Figure 
26b shows the 500 samples of the waveform in (a) after filtering with 
a nonrecursive filter of the type shown in Fig. 24, whose gain was zero 
after 3.2 kHz. Figure 26c shows 160 samples of the result of a change 
of sampling rate from 20 kHz to 6.4 kHz. The value of N was 700 and 
N' = (6,400) (700)/20,000 = 224. It is difficult to judge quantitatively 
from such a figure the accuracy of the interpolation. It does seem safe 
to conclude that the error is not extreme. Our experience has been that 
there is significant error only in the first and last few samples of the 
N' output samples. Using the chirp z-transform algorithm, these "bad" 
samples need not ever be computed, for example, only M "good" values 
need be computed. 




10 15 

TIME IN MILLISECONDS 



Fig. 26 — An example of interpolation for the purpose of changing the sampling 
rate: (a) 500 samples of speech at 20 kHz sampling rate; (b) 500 samples of (a) 
after low pass filtering to 3.2 kHz; (c) 160 samples of (a) after changing the 
sampling rate to 6.4 kHz using the chirp z-transform. (In all cases the samples 
are connected by straight lines). 



CHIRP Z-TRANSFORM 1287 

If one wishes to low-pass filter a waveform and then go to a lower 
sampling rate, the filtering and interpolation can be combined if we 
use a nonrecursive filter. That is, the discrete Fourier transform of the 
filter impulse response can be simply combined with H N (k). 

V. LIMITATIONS 

Several times we have pointed out shortcomings of the chirp z-trans- 
form algorithm. One limitation in using it to evaluate the z-transform 
off the unit circle stems from the fact that we may be required to com- 
pute Wt n ' /2 for large n. If W differs very much from 1.0, TP*" V2 can 
become very large or very small when n becomes large. (We 
require a large n when either M or N become large, since we 
need to evaluate W*' /a for n in the range — N < n < M.) For example, if 
W = e - 2 - 5/l0 - 000 « 0.999749, and n = 1,000, JP*" va = e il28 which 
exceeds the single precision floating point capability of most computers 
by a large amount. Hence the tails of the functions W ±n ' /2 can be 
greatly in error, thus causing the tails of the convolution (the high fre- 
quency terms) to be grossly inaccurate. The low frequency terms of the 
convolution will also be slightly in error, but these errors generally are 
negligible. 

An example of this effect is shown in Fig. 27. The contour for the 
five curves in this figure was held fixed (contour 5 in Fig. 6) and the 
number of frequency points in the range to 5,000 Hz was increased 
in steps of 2 from 50 to 800. Spectral samples are plotted every 100 Hz 
for comparison. (This example was programmed using single-precision 
floating-point arithmetic on a GE 635 computer with a 36 bit word 
length.) It is seen that as the number of output points increases, errors 
in the high frequency region become large and completely mask the 
fifth resonance for the 800 point case. The effects of the inaccuracy in 
W ±n ' /2 can also be seen at low frequencies. For example, the spectral 
magnitude at Hz goes from about 120 dB to 134 dB as the number of 
points goes from 50 to 800. These small errors generally do not affect the 
gross spectral characteristics as seen in Fig. 27. The resonances are 
easy to locate in all cases until the errors get exceedingly large. One 
can push the maximum point limit higher than 800 (in this case) by 
using double precision arithmetic. 

The limitation of contour distance in or out from the unit circle is 
again the result of computation of W ±n ' /2 . As W deviates significantly 
from 1.0, the number of points for which W ±n ' /2 can be accurately 
computed decreases. It is of importance to stress, however, that for 



1288 



THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 




1000 2000 3000 4000 
FREQUENCY IN HZ 



Fig. 27 — A comparison of magnitude plots for varying number of points on the 
same spiral contour. The fifth plot shows the effect of errors in evaluating W n2/J for 
large n. (Points are plotted every 100 Hz in each curve to aid in comparison.) 

W = 1 there is no limitation of this type since W ±n ' /2 is always of 
magnitude 1. 

The other main limitation of the chirp z-transform algorithm stems 
from the fact that two L point fast Fourier transforms and one L/2 
point fast Fourier transform must be evaluated where L is the smallest 
convenient integer greater than N + M — 1 as previously mentioned. 
We need one fast Fourier transform and 2L storage locations for the 
transform of x n A~ n W n ' /2 ; one fast Fourier transform and L+2 storage 
locations for the transform of W~ n ' /2 ; and one fast Fourier transform 
for the inverse transform of the product of these two transforms. We 
do not know a way of computing the transform of W~ n ' /a either re- 
cursively or by a specific formula (except in some trivial cases.) Thus 
we must compute this transform and store it in an extra L + 2 storage 
locations. Of course, if many transforms are to be done with the same 
value of L we need not compute the transform of W~ n ' /2 each time. 

We can compute the quantities A~ n W n7/2 recursively, as they are 



CHIRP Z -TRANSFORM 1289 

needed, to save computation and storage. This is easily seen from the 
fact that 

A- b+ V» ,fl - (A-W^-WlflA- 1 . (45) 

If we define 



and 



then 



and 



C n = A-W n ' /2 (46) 

D n - WW+A- 1 (47) 

D n+l = W-D n (48) 

C n+1 = C n -D„. (49) 

Setting A = 1 in equations (45) through (49) provides an algorithm 
for the coefficients required for the output sequence. A similar recursion 
formula can be obtained for generating the sequence A~ n W in 
The user is cautioned that recursive computation of these coefficients 
may be a major source of numerical error, especially when W ~ 1 
or <p tt 0. 

VI. SUMMARY 

We give a computational algorithm for numerically evaluating the 
2-transform of a sequence of A r time samples. This algorithm, the 
chirp 2-transform algorithm, enables the evaluation of the z-trans- 
form at M equiangularly spaced points on contours which spiral in or 
out (circles being a special case) from an arbitrary starting point in 
the z-plane. In the s-plane the equivalent contour is an arbitrary 
straight line. 

The chirp 2-transform algorithm has great flexibility in that neither 
A r or M need be composite numbers; the output point spacing is 
arbitrary ; the contour is fairly general and N need not be the same as 
M. The flexibility of the chirp 2-transform algorithm comes from being 
able to express the 2-transform on the above contours as a convolu- 
tion, permitting the use of well-known high speed convolution tech- 
niques to evaluate the convolution. 

Applications of the chirp 2-transform algorithm include enhance- 
ment of poles for use in spectral analysis, high resolution narrowband 



1290 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

frequency analysis, and time interpolation of data from one sampling 
rate to any other sampling rate. These applications are explained in 
detail. The chirp 2-transform algorithm also permits use of a radix 2 
fast Fourier transform program or device to compute the discrete 
Fourier transform of an arbitraiy number of samples. Examples were 
presented illustrating how the chirp z-transform algorithm was used 
in specific cases. It is anticipated that other applications will be found. 

VII. ACKNOWLEDGMENT 

The authors would like to acknowledge the help and assistance of 
Miss Carol McGonegal in programming many of the examples shown 
in this paper and Miss Debbie Hougak for valuable clerical help in 
preparing the manuscript. 

APPENDIX 

Fast Fourier Transforms for Two Real L Point Sequences 

The purpose of this appendix is to show how the fast Fourier trans- 
forms of two real, symmetric L point sequences can be obtained using 
one L/2 point fast Fourier transform. 

Let x„ and y n be two real, symmetric L point sequences with corre- 
sponding discrete Fourier transforms X k and Y k . By definition, 

x n = x L - n n = 0,1,2, ■•• ,L- 1; 

Vn = Vt-n 

it is easily shown that X k and Y k are real, symmetric L point se- 
quences, so that 

X k = X h . k k = 0lt2t ... ,L- 1. 
Y k = Y L - t 

Define a complex, L/2 point sequence Un whose real and imaginary 
parts are 



Re [U„] = X 2n - ?/2n + l + V2n-1 

Im [u„] = y 2n + x 2n+1 — x 2n -i. 



n = 0, 1, ••• ,1/2-1. 



The L/2 point discrete Fourier transform of u n is denoted U k and is 
calculated by the fast Fourier transform. The values of X k and Y k may 



CHIRP Z-TRANSFORM 1291 

be computed from U k using the relations 
X k = %{Re[U k ] + Re[U L/2 . k ]\ 



1 



a 2ir , 

4 sin -=- k 

Li 



[Re[U k ] - Re[U L/2 . k }\ 



Y k = i{Tm[U k ] + lm[U W2 . k ]\ 

1 



4 sin -=- k 

Li 



\lm[U k \ - lm[U L/2 . k )} 



for k = 1,2, ••■ ,L/2 - 1. 
The remaining values of X k and Y k are obtained from the relations 



n-0 
L-l 

Y =TiVn 

n-0 

x L/2 = S (- dx 

n-0 



REFERENCES 



1. Cooley, J. W. and Tukey, J. W., "An Algoritiim for the Machine Calculation 

of Complex Fourier Series," Mathematics of Computation, 19, No. 90 
(April 1965), pp. 297-301. 

2. "What is the Fast Fourier Transform?," G-AE Subcommittee on Measure- 

ment Concepts, IEEE Tran. Audio and Electroacoustics, AU-16, No. 2, 
(June 1967), pp. 45-55. 

3. Stockham, T. G., "High Speed Convolution and Correlation," 1966 Spring 

Joint Computer Conference, Amer. Federation of Inform. Processing Soc. 
Proc, 28, Washington, D. C, April 1966, pp. 229-233. 

4. Helms, H. D. "Fast Fourier Transform Method of Computing Difference 

Equations and Simulating Filters," IEEE Trans. Audio and Electroacous- 
tics, AU-16, No. 2 (June 1967) , pp. 85-90. 

5. Bluestein, L. I., "A Linear Filtering Approach to the Computation of the 

Discrete Fourier Transform," 1968 Northeast Electronics Research and 
Engineering Meeting Record, 10, November 1968, pp. 218-219. 

6. Rader, C. M. and Gold, B., "Digital Filter Design Techniques in the Fre- 

quency Domain," Proc. IEEE, 55, No. 2, (February 1967), pp. 149-171. 



1292 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969 

7. Schafer, R. W. and Rabiner, L. R., Automatic Formant Analysis of Speech 

Using the Chirp z-Transform Algorithm" to be presented at the 1969 IEEE 
International Conference on Communications, Boulder, Colorado, June, 
1969. 

8. Gold, B. and Jordan, K. L., "A Direct Search Procedure for Designing Finite 

Duration Impulse Response Filters," IEEE Trans, on Audio Electro- 
acoustics, AU-17, No. 1 (March 1969) pp. 33-36. 

9. Gentleman, W. M. and Sande, G., "Fast Fourier Transforms for Fun and 

Profit," 1966 Fall Joint Computer Conference, American Federation of 
Information Processing Societies Proc, 29, Washington, D. C, November 
1966, pp. 563-578.